program example_MySparseOperator
	use MySparseOperator
	implicit none
	real(kind=8), dimension(:), allocatable :: a
	integer, dimension(:), allocatable :: ia, ja
	real(kind=8), dimension(:), allocatable :: b, x
	integer, dimension(:), allocatable :: ib,jb

	integer :: mtype, msglvl



	allocate(a(17)) 
	a = (/&
		&	1.d0,	2.d0,	&
		&		3.d0,	4.d0,	&
		&			5.d0,	6.d0,	&
		&				7.d0,	8.d0,	&
		&					9.d0,	10.d0,	&
		&						11.d0,	12.d0,	&
		&							13.d0,	14.d0,	&
		&								15.d0,	16.d0,	&
		&									17.d0	/)
	allocate(ia(10))
	ia = (/ 1,3,5,7,9,11,13,15,17,18 /)
	allocate(ja(17))
	ja = (/&
		&	1,	2,	&
		&		2,	3,	&
		&			3,	4,	&
		&				4,	5,	&
		&					5,	6,	&
		&						6,	7,	&
		&							7,	8,	&
		&								8,	9,	&
		&									9	/)
	write(*,*) '============== sparse matrix A to be deleted rows or columns============='
	call mso_prtSpMtxR(a, ia, ja, (/1,9/), (/1,9/), 'f10.1')

	! 	delete some rows of a CSR format sparse matrix
	call mso_dcsrdltRC(a,ia,ja,(/2,2,5/),'R', b,ib,jb)
	write(*,*) '============== sparse matrix B from A which rows are deleted ============='
	call mso_prtSpMtxR(b, ib, jb, (/1,size(ib)-1/), (/1,maxval(jb)/), 'f10.1')


	if ( allocated(a) ) deallocate(a)
	if ( allocated(ia) ) deallocate(ia)
	if ( allocated(ja) ) deallocate(ja)  
	allocate(ia(8 + 1))
	ia = (/ 1, 5, 8, 10, 12, 15, 17, 18, 19 /)
	allocate(ja(18))
	ja = (/ 1,    2,       6, 7,    &
           2, 3,    5,          &
              3,             8, &
                 4,       7,    &
                    5, 6, 7,    &
                       6,    8, &
                          7,    &
                             8 /)
	allocate(a(18))
	a = (/ 7.d0,        1.d0,             2.d0, 7.d0,        &
             -4.d0, 8.d0,       2.d0,                    &
                    1.d0,                         5.d0,  &
                          7.d0,             9.d0,        &
                                5.d0, 1.d0, 5.d0,        &
                                     -1.d0,       5.d0,  &
                                           11.d0,        &
                                                  5.d0 /)
	allocate(x(8))
	x = 1.d0	! expected solution
	if ( allocated(b) ) deallocate(b) 
	allocate(b(8))
	call mkl_dcsrsymv('u', 8, a, ia, ja, x, b)
	mtype = -2	! symmetric indefinite
	msglvl = 1	! show all message data
	call mso_dpardiso(mtype,msglvl,a,ia,ja,b,x)


	if ( allocated(a) ) deallocate(a)
	if ( allocated(ia) ) deallocate(ia)
	if ( allocated(ja) ) deallocate(ja)
	if ( allocated(b) ) deallocate(b)
	if ( allocated(ib) ) deallocate(ib)
	if ( allocated(jb) ) deallocate(jb)
	if ( allocated(x) ) deallocate(x) 
end program example_MySparseOperator